On the saturation amplitude of the f-mode instabihty 
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We investigate strong nonlinear damping eflects which occur during high amplitude oscillations 
of neutron stars, and the gravitational waves they produce. For this, we use a general relativistic 
nonlinear hydrodynamics code in conjunction with a fixed spacetime (Cowling approximation) and a 
polytropic equation of state (EOS) . Gravitational waves are estimated using the quadrupole formula. 
Our main interest are I — m = 2 f modes subject to the CFS (Chandrasekhar, Friedman, Schutz) 
instability, but we also investigate axisymmetric and quasiradial modes. We study various models 
to determine the influence of rotation rate and EOS. We find that axisymmetric oscillations at high 
amplitudes are predominantly damped by shock formation, while the nonaxisymmetric / modes are 
mainly damped by wave breaking and, for rapidly rotating models, coupling to nonaxisymmetric 
inertial modes. From the observed nonlinear damping, we derive upper limits for the saturation 
amplitude of CFS-unstable / modes. Finally, we estimate that the corresponding gravitational 
waves for an oscillation amplitude at the upper limit should be detectable with the advanced LIGO 
and VIRGO interferometers at distances above 10 Mpc. This strongly depends on the stellar model, 
in particular on the mode frequency. 



PACS numbers: 04.30.Db, 04.40.Dg, 95.30.Sf, 97.10.Sj 

I. INTRODUCTION 

As a consequence of Einstein's theory of relativity, it 
is predicted that the violent nonradial pulsations excited 
immediately after the birth of a rotating proto-neutron 
star result in the emission of significant amounts of grav- 
itational radiation. In addition, it has been discovered 
by Chandrasekhar [T and Friedman & Schutz [51 13] that 
certain pulsation modes of rotating relativistic stars may 
grow exponentially due to gravitational radiation back- 
reaction, especially if the proto-neutron star is rapidly 
rotating; this is the so-called CFS instability. 

The exact amount of energy emitted by the stable os- 
cillations of neutron stars in the form of gravitational 
waves is actually unknown and depends uniquely on the 
initial conditions. The stable stellar pulsations will emit 
gravitational waves which are detectable only if the neu- 
tron star is in our galaxy or in the most optimal case 
in the nearby ones [4]. Thus, the cases in which these 
oscillations become unstable are of special interest for 
gravitational wave astronomy, since the corresponding 
gravitational waves will be detectable even for sources 
in the nearby galactic clusters. It has been suggested 
that the detection of gravitational waves from pulsating 
neutron stars will allow the study of their interior, see 
[SH?]. It is expected that the identification of specific 
pulsation frequencies in the observational data will re- 
veal the properties of matter at densities that currently 



cannot be probed by any other experiment. 

The study of the dynamics of fast rotating neutron 
stars in a general relativistic framework was practically 
impossible until recently. The main reason why linear 
theory failed is the form and the size of the relevant per- 
turbation equations. Thus, it is not surprising that the 
first results for the oscillations of fast rotating stars were 
derived by using evolutions of the nonlinear equations, 
see [5HTT|. Still most of these studies were purely ax- 
isymmetric and thus there was a significant influence of 
rotation on the oscillation spectra only for very high ro- 
tation rates. The CFS instability on the other hand is ac- 
tive only for nonaxisymmetric perturbations. In the last 
two years there was a significant progress in the study 
of nonaxisymmetric perturbations, for the first time in a 
general relativistic framework, using both perturbation 
theory [71 [T^HTi] and nonlinear evolutions of coupled hy- 
drodynamic and Einstein equations [15] . 

For the case of r modes, the CFS instability is active 
for any rotation rate [T5H20) . In practice the unstable 
r modes will be excited even for slowly rotating stars if 
they are hot enough (T > 10^ K), so that the instability 
growth time is shorter than the shear viscosity damping 
time. Initially it was considered as a prime source for 
gravitational waves, but later more detailed studies on 
the effect of the magnetic field on the instability [21], or 
the presence of hyperons in the core [2^ , seriously ques- 
tioned the potential of the instability. Detailed stud- 
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ies suggested that the r mode is limited to very small 
amplitudes due to energy transfer to a large number of 
other inertial modes, in the form of a cascade, leading 
to an equilibrium distribution of mode amplitudes |23j . 
The small saturation values for the amplitude are sup- 
ported by recent nonlinear estimations |24j |25] based on 
the drift, induced by the r modes, causing differential ro- 
tation. On the other hand, hydrodynamical simulations 
of limited resolution showed that an r mode of large ini- 
tial amplitude does not decay appreciably over several 
dynamical timescales [26j . but on a longer timescale a 
catastrophic decay occurs [17] ; this indicates a transfer of 
energy to other modes due to nonlinear mode couplings, 
and suggests that a hydrodynamical instability may be 
operating. A specific resonant 3-mode coupling was iden- 
tified in [28i as the cause of the instability and a pertur- 
bative analysis of the decay rate suggests a maximum 
dimensionless saturation amplitude Ofmax < 10~^ — 10"'^. 
A new computation using second-order perturbation the- 
ory finds that the catastrophic decay seen in the hydro- 
dynamical simulations |27l I28j can indeed be explained 
by a parametric instability operating in 3-mode couplings 
between the r mode and two other inertial modes . 

The CFS instability of the / modes (nonradial modes 
with no radial nodes, restored by pressure) can poten- 
tially be the strongest source for gravitational waves from 
isolated neutron stars. Although it was the first to be 
studied, already 30-years ago, it was forgotten because 
at that time it was nearly impossible to study rapidly ro- 
tating neutron stars in a general relativistic framework. 
Studies using Newtonian theory came to the disappoint- 
ing conclusion that (in Newtonian theory) the I = m = 2 
f mode is stable for rotation rates below the Kepler limit. 
The I = m = A f mode does become unstable, but the 
growth time is so long that viscosity severely limits the 
amplitudes, rendering the process unimportant for grav- 
itational wave astronomy. Only in the late '90s, Ster- 
gioulas & Friedman [34] showed that for certain EOS the 
I = ni — 2 f mode can become unstable in full GR (gen- 
eral relativity). Still, this study was left aside since no 
proper linear or nonlinear codes were available to deal 
with the perturbations of rapidly rotating neutron stars 
in GR. 

If such an instability sets in, the star will emit copious 
amounts of gravitational radiation. First studies suggest 
that the signal may be detectable from distances as far as 
15 Mpc for the current sensitivity of Virgo & LIGO, and 
from distances greater than 100 Mpc for the sensitivities 
of Advanced Virgo & LIGOs [55H37) . In the latter case 
the event rate of supernovas resulting in the creation of 
a proto-neutron star can be quite high, i.e. more than 
thousand per year. However, it should be noted that it 
is not clear what will be the initial rotation period of the 
collapsing core; it seems to depend strongly on the profile 
of the angular momentum distribution, the initial mass 
and angular momentum of the collapsing star, and also 
the strength of the magnetic field, which can slow down 
the newly born compact object quite efficiently. The 



hope is that still a number of supernova events, maybe 
10 %, will produce rapidly rotating cores which can be 
subject to the CFS instability of the / modes. In any 
case, the next generation of gravitational wave detectors, 
such as Einstein Telescope (FT) ^38j will be ideal for the 
detection of this type of instability [4j . 

The two most recent hydrodynamical simulations [361 
I37j (in the Newtonian limit and using an artificially in- 
creased post-Newtonian radiation-reaction potential) es- 
sentially confirm this picture. In [36' a differentially 
rotating, N = 1 polytropic model with a large /? = 
r/|VK| ~ 0.2... 0.26 was chosen as the initial equilib- 
rium state (r is the rotational kinetic energy and \W\ 
the gravitational binding energy). The main difference 
of this simulation compared to the ellipsoidal approxima- 
tion [3S] lies in the choice of the EOS. For N = 1 Newto- 
nian polytropes it was argued that the secular evolution 
cannot lead to a stationary Dedekind-like state. Instead, 
the /-mode instability will continue to be active until all 
nonaxisymmetries are radiated away and an axisymmet- 
ric shape is reached. In another recent simulation J37j, 
the initial state was chosen to be a uniformly rotating, 
N = 0.5 polytropic model with T/\W\ ~ 0.18. Again, 
the main conclusions reached in [35' were confirmed, how- 
ever the assumption of uniform initial rotation limits the 
available angular momentum that can be radiated away, 
leading to a detectable signal only out to about 40 Mpc. 
The star appears to be driven towards a Dedekind-like 
state, but after about 10 dynamical periods, the shape 
is disrupted by growing short-wavelength motions, which 
are explained by a shearing type instability such as the 
elliptic flow instability [5^ . 

The recent progress [71 [T2l - fT5| , demonstrates that the 
problem can be handled properly, which should allow to 
answer the following questions: First, has the unstable 
mode any chance to be excited? Second, what is the 
instability window of the mode? Third, what are the 
exact rotation rates at the onset of the instability for 
the various EOS? And finally, what is the maximum am- 
plitude that the mode may reach before saturated by 
nonlinear phenomena, e.g. mode coupling, mass loss, or 
shock formation? In this work we try to answer the last 
question. In particular, we investigate nonlinear effects 
active at oscillation amplitudes sufficiently large to allow 
a detection of sources at distances greater than 1 Mpc. 
Since our methods are not restricted to a certain mode, 
we map out nonlinear effects occurring for axisymmetric 
modes as well. The latter also serves as a numerically 
less expensive test case to calibrate our code. 

The question about the saturation amplitude of the 
/-mode instability is a major one that has to be an- 
swered since it may completely eliminate the importance 
of this instability. In this article we present a systematic 
study of nonlinear damping effects for different ampli- 
tudes, rotation rates, and equations of state. For this, we 
use a nonlinear general relativistic evolution code named 
PIZZA [ini [H] . Up to now it has been successfully tested 
on small amplitude stable oscillations of rapidly rotating 
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neutron stars and self-gravitating tori. After a number of 
improvements the code is well suited to follow even large 
amplitude oscillations of rotating stars, such as those ex- 
cited during the saturation phase of the /-mode instabil- 
ity. A similar study has been performed for the r-mode 
instability in [27], where an artificially excited r mode 
was left to evolve. The amplitude of the r mode seemed 
initially to decrease slowly but then it decayed catas- 
trophically. The decay time was found to be amplitude 
dependent and the effect attributed to three mode cou- 
pling. 

The structure of the paper is as follows. First we de- 
scribe the numerical methods in Sec. ITTl In Sec. 
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introduce the stellar models and physical problems we 
studied. The numerical results are presented in Sec. |IV[ 
The main results on the damping and detectability of 



the CFS-unstable / mode can be found in Sec. IV E and 
IV H and a summary is given in Sec. [Vj 

Throughout the paper we use geometrical units c — 
G = 1. Greek indices run from ... 3, Latin indices from 
1 . . . 3. We denote the lowest order quasi radial mode 
by -F, the axisymmetric fundamental pressure modes 
(/ modes) by %, and nonaxisymmetric ones by '/m, 
where to > refers to the counter-rotating modes. 



II. NUMERICAL METHOD 
A. Time evolution 

For our simulations, we use the pizza code described 
in [To]. It evolves the general relativistic nonlinear hy- 
drodynamic evolution equations for an ideal fluid without 
magnetic fields, while keeping the spacetime metric fixed 
(Cowling approximation) . 

The hydrodynamic equations in covariant form are 



0, 



(pu^) = 0. 



(1) 
(2) 



The stress energy tensor T^'^ of an ideal fluid is given by 



T^"^ = phu'^u" + Pg- 



(3) 



where p, ft., and are the rest mass density, relativistic 
specific enthalpy, and 4-velocity of the fiuid; g^'^ is the 
metric tensor of signature (—,+,+,+). 

The covariant equations can be written as a first order 
system of hyperbolic evolution equations in conservation 
form with source terms: 



q={D,E,Sj). 



(4) 
(5) 



For our study, it is important to use a conservative formu- 
lation, because it ensures correct shock wave propagation 
speeds when evolved by means of a finite volume scheme. 



The evolved hydrodynamic variables are 
D = ^Wp, 

E = ^ {W^ph -P- Wp) , 
= ^W'^phv,, 



(6) 
(7) 
(8) 



where 7 is the determinant of the 3-metric, the 3- 
velocity, and W the corresponding Lorentz factor. In the 
Newtonian limit using Cartesian coordinates, D, E, and 
reduce to mass density, energy density, and momen- 
tum density. 

The flux terms P = {fhJh^ fs,) are given by 



fh = ^'D, 

fh = w'E + a^v'P, 
Ps^=w^S,+a^P5], 



(9) 
(10) 
(11) 



where = /vP is the coordinate velocity of the fluid 
and a the lapse function. The evolution equations need 
to be completed by an equation of state (EOS) to com- 
pute the pressure. We choose a polytropic EOS defined 

by 



P[p) = Pp 



Pp 



Kp\ 



(12) 



where pp is a constant density scale, and the constant F 
is the polytropic exponent. 

When assuming a one-parametric EOS, the evolution 
equations ^ are over-determined. Therefore, we do not 
evolve the energy density E, but recompute it from the 
remaining variables and the EOS. The physical implica- 
tions will be discussed in the next subsection. 

The above formulation is used by many modern rela- 
tivistic codes. For a review, we refer to [4(1. Like most 
codes, the pizza code is based on a HRSC (high reso- 
lution shock capturing) scheme, which was however op- 
timized for quasistationary simulations. The difference 
to standard HRSC schemes is that source and fiux terms 
are treated more consistently, making use of a special 
formulation of the source terms in Eq. Q. As a conse- 
quence, the code is able to preserve a stationary star to 
high accuracy. For details, see [10] . 

One of the weak points of current relativistic hydro- 
dynamic schemes is the treatment of the stellar surface, 
where the numerical schemes designed for the interior 
would break down without further measures. The most 
common remedy is to enforce an artificial atmosphere of 
low density. The original pizza code described in [TO] 
used a different workaround, which works perfectly for 
low amplitude oscillations. Unfortunately, we found that 
it yields unsatisfactory results for high amplitude oscil- 
lations in conjunction with a stiff EOS. 

For the results in this article, we used yet another 
method of treating the surface. In short, it is based on a 
smooth transition of the numerical fiux, between the for- 
mally correct expression at and above a certain density. 
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and a flux which causes no acceleration at zero density. 
The local density scale for the beginning of the transi- 
tion is computed from the local gravity, using an expres- 
sion which corresponds to the average density gradient 
of the stationary model near the stellar surface times 
the grid spacing. No artificial atmosphere is used. To 
prevent negative densities, the flux is limited. The de- 
tails of the scheme will be described in a forthcoming 
publication. We just note that the evolution below the 
aforementioned density scale is still plain wrong, like it 
is the case with other common schemes. Due to the low 
densities involved, the huge local errors at the surface 
have only limited impact on the global evolution. Never- 
theless, the surface is an important source of numerical 
errors for this study, in particular with regard to the nu- 
merical damping. 

The code has been tested on shock tube problems, lin- 
ear oscillations of neutron stars, see [101 HH SI]i and 
recently with self-gravitating relativistic tori. On the 
timescales of the simulations presented here, the code 
is able to evolve the stationary background model with- 
out significant changes in density or rotation profile. In 
contrast to codes using artificial atmospheres, our code 
exactly conserves the total mass. 

We use Cartesian grids for three-dimensional simula- 
tions, cylindrical coordinates for two-dimensional ones, 
and spherical coordinates for one-dimensional problems. 
The rigidly rotating models are evolved in the corotat- 
ing frame. To save computation time, we further assume 
equatorial symmetry. Finally, we apply vacuum outer 
boundary conditions, such that any material hitting the 
outer boundary is lost, and monitor the total mass to 
detect this case. 



B. Treatment of shock waves 

For some of our results, shock formation plays an im- 
portant role. The use of a polytropic EOS then becomes 
one of the main limitations, as will be explained in the 
following. 

For any ideal fluid system with a two-parametric EOS, 
the physically correct solution evolves adiabatically as 
long as there are no shock waves present. For initial 
data which is isentropic, the evolution in the absence of 
shocks is thus the same as for an one-parametric EOS 
corresponding to a curve of constant specific entropy of 
the two-par amctric EOS. 

The polytropic EOS used in our simulations are the 
curves of constant specific entropy of the ideal gas EOS 
given by 



P(p,e) = (r-l)pe, 



(13) 



where e is the specific internal energy and F matches 
the polytropic exponent. Our simulations are therefore 
equivalent to the case of an ideal gas EOS with isentropic 
initial data, until shock waves form. Hence we can accu- 
rately predict the onset of shock formation. 



After the shocks form, using the polytropic EOS be- 
comes unphysical. Formally, the reduced system of evo- 
lution equations for the polytropic case admits solutions 
with discontinuities. However those differ from realistic 
shock waves in two aspects: there is no entropy produc- 
tion, i.e. shock heating, and the local conservation of 
energy is violated. 

Nevertheless, we also present results based on the evo- 
lution after shock formation. For this, we make the as- 
sumption that the main difference between shock solu- 
tions for hot and cold EOS is that the kinetic energy 
converted into heat for the correct solution is simply lost 
when using the cold EOS, while the evolution of the den- 
sity profile remains similar. Although this is clearly a 
strong simplification, the results may serve as a first es- 
timate. 

The reason why we do not simply use the ideal gas 
EOS is purely technical: the version of our code which 
is capable of evolving two-parametric EOSs is using a 
method of treating the surface which is not well suited 
for high amplitude oscillations. Fortunately, the only re- 
sults depending on evolutions beyond shock formation 
are the damping times of axisymmetric oscillations, but 
not the critical amplitudes for the onset of shock forma- 
tion; for the nonaxisymmetric oscillations, no significant 
shocks occur, and hence the use of the polytropic EOS is 
perfectly valid there. 



C. Estimating the gravitational wave amplitude 

To extract gravitational waves, we use the multipole 
formalism for Newtonian sources, as described in (42] . 
In this formalism, the radiation field at distance r far 
away from the source is expanded in terms of spin tensor 
harmonics 



-nE2.ln 



1=2 m=-l 

and the total gravitational luminosity is given by 

oo I 



Ini I 



Im I 



(14) 



(15) 



1=2 m=-l 



In the Newtonian limit, the coefficients of the radiation 
field multipole expansion can be expressed in terms of 
the mass- and current-multipoles of the source 



aE2 _ 
^Irn — 



aB2 _ 
^Im — 



4x/2^ {l + l){l + 2) 
(2Z + 1)!!V (1-1)1 

32n 



(21 + I) dlqim, (16) 



{21 + iyAy 2(/-i) 



(17) 
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with the multipole moments given by 



2/ + 1 
■hm = I pr'-v-Yg 



K-t vim* ^3^^ 



(18) 
(19) 



The above integrals are defined in the incrtial frame, and 
Y^"' denotes the vector spherical harmonics of magnetic 
type. 

The error induced by using the above formulas for 
sources as relativistic as neutron stars is difhcult to esti- 
mate analytically. In [33], results using the quadrupole 
formula are compared with direct wave extraction meth- 
ods for the case of a pulsating neutron star model, finding 
an error around 10-20 % in the strain amplitude. 

However, this error does not directly carry over to our 
results. The reason is that besides neglecting relativistic 
corrections, another error arises from the ambiguity of 
choosing a coordinate system for strongly curved space- 
times. This depends not only on the model, but also on 
the gauge choices made when computing the initial data. 
To get a rough estimate for the deviation from Euclidean 
geometry, we compute 



m = 



V3 



■pe 



R 

Re 
Rci 



R„ 



VP 



R 

Rp 



(20) 
(21) 



where Re and Rp are the equatorial and polar coordinate 
radius, while Rpe, Rpp, and i?circ are the proper equato- 
rial radius, proper polar radius, and equatorial circum- 
ferential radius. 

From this average quantities we estimate the error of 
the multipole moments by assuming that the radius r in 
Eq. ( 18 1 and Eq. ( 19 ) is wrong by a constant factor rjr, in 
the sense that the correct value is in the range [r/rjr, rr]r], 
and that the volume element d^x is wrong by a factor of 
7^4. We further set r]r to the maximum over 771... 3 and the 
reciprocal values. 

In consequence, the strain amplitude (which is mainly 
due to the / = 2 mass multipoles) would be wrong by 
a factor rjt = 77^ '74- While for most of our models, rit is 
in the range 1.2-2, it is as big as 13 in one case. This 
shows how dangerous it is to generalize error bounds for 
the quadrupole formula computed for a "typical" neutron 
star model. 

Our estimate does not take into account the possibility 
that large cancellations occur in the integrals for the mul- 
tipole moments. In this study, this affects the quasira- 
dial F-mode oscillations, in particular for slowly rotating 
stars. For those, we do not compute the error; to obtain 
robust results, fully relativistic studies are needed. 

To implement the above formalism, we compute the 
multipole moments in the corotating frame used in our 
simulations. For the current multipoles, we neverthe- 
less use the incrtial frame velocity in the integrals. It is 



straightforward to show that the multipole moments in 
the incrtial frame are given by 



J —imilt jf 

'Jim — 'J] 



(22) 
(23) 



where and are the multipole moments computed 
in the corotating frame defined hy (j)' = (j) — fit, with 
being the angular velocity of the star. 

The multipole moments are computed during the evo- 
lution with a sampling rate of 75 kHz. When evaluating 
time derivatives of 2nd order or higher, special care has to 
be taken to avoid amplification of high frequency numeri- 
cal noise. For this, we first apply numerical smoothing by 
convolution of the time series with a Blackman window 
function. Derivatives up to 3rd order are then computed 
using cubic splines. For derivatives of 4th and 5th order, 
we first compute the 3rd derivative and than apply the 
whole scheme again. 

We computed the frequency response function of the 
resulting scheme, which effectively cuts off contributions 
to the gravitational wave (GW) signal above 10 kHz. In 
the frequency range of the actual signal, the loss of strain 
amplitude due to the smoothing stays below 15 %. 

Although only the I — 2 mass multipoles are important 
for our problems, we generally compute all mass multi- 
pole moments up io I — A, and the current multipoles 
with m > 2, / < 4, 

The second important error for the strain, beside the 
use of the multipole formula, is due to the Cowling ap- 
proximation, which is known to cause a significant error 
in the oscillation frequencies. Given an estimate for the 
incrtial frame frequency fi of a harmonic oscillation in 
full GR, one can approximate the strain amplitude 
by 



A 



E2 



A. 



E2 



(24) 



where fi and are frequency and strain amplitude in 
the Cowling approximation. 



D. Measuring dissipation 

In order to measure the damping of oscillations, we 
monitor the evolution of an average corotating velocity 
defined by 



(25) 
(26) 
(27) 
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where a and /3* are lapse function and shift vector. Since 
we are working in the corotating frame, is a corotating 
velocity. 

This measure is zero if and only if the fluid velocity is 
everywhere the same as for the unperturbed stationary 
star. For the nonrotating case in the Newtonian limit, v 
is also directly related to the total kinetic energy. As long 
as there is only one dominant oscillation mode, the decay 
of is a measure for the decay of the mode amplitude. 
On the other hand, it is not sensitive to energy transfer 
from one oscillation mode to another. 

From the evolution, we extract three quantities. 

• The initial amplitude Aj = v(0). 

• The final amplitude. Since v is oscillating, we use 
the maximum amplitude during the last oscillation 
cycle. 



Af = max{v{t)\T, -To<t< T,}, 



(28) 



where Tf, is the time over which we evolved the 
system and Tq is the period of the oscillation mode 
used to perturb the system. 

• The timescale r of the initial decay of w, which we 
define as 



VaiO) ' 



ft+To 



dt, (29) 



where Tq is the oscillation period of the mode used 
for perturbation, and AT — ATq. 

To detect the presence of shocks, we make use of the 
fact that there exists a conserved energy when using the 
Cowling approximation, given by 



Pe = -ar° - Wp. 



(30) 
(31) 



More precisely, Ec is conserved for smooth and weak so- 
lutions of Eq. ([T]) . As mentioned, shock solutions of the 
reduced system of equations used with a one-parametric 
EOS violate the local energy conservation, and thus the 
conservation of Ec. Thus, any change of E^ points to the 
existence of shock waves. 



E. Computing eigenfunctions and eigenfrequencies 

To extract eigenfunctions, we use the mode recycling 
method described in [9|. In short, we perturb the star us- 
ing a generic perturbation to excite different modes and 
extract their frequencies using Fourier analysis. Then we 
extract a first guess of the eigenfunction by evaluating 
at each point of the numerical grid the Fourier integral 
of specific energy and velocity perturbations at the fre- 
quency of the desired mode. 



The estimate of the eigenfunction obtained in this way 
is in general still contaminated with other modes, due to 
the finite evolution time. Therefore, additional simula- 
tions are performed, using the eigenfunction obtained in 
the previous step as initial perturbation, until only the 
desired oscillation mode is present in the evolution. 

We further improved this scheme by making use of the 
fact that for any axisymmetric star, the eigenfunctions 
of e,v'^, w^, and can be written in the form 



(32) 



The two-dimensional eigenfunction Si{d, z) is real-valued 
and has a well defined z-parity. The constant /3 is an 
irrelevant complex phase. 

The Fourier analysis of the numerical evolution yields 
an estimate for the complex-valued bt(d^ z, 0). To obtain 
(5e, we first compute the integral 



-im(f} 



(33) 



Since we use Cartesian coordinates for nonaxisymmet- 
ric problems, this step involves interpolation. For this, 
a multidimensional cubic spline interpolation method is 
applied. Note we also increase the resolution by evaluat- 
ing the above integral at finer intervals than the original 
grid spacing. 

Errors due to the presence of unwanted modes with dif- 
ferent (/)-dependency are strongly suppressed during the 
computation of Sic- This step is absolutely necessary to 
separate co- and counter-rotating modes for slowly rotat- 
ing stars, since their frequencies approach each other. 

It seems plausible to take the real or imaginary part 
of Sec as estimate for the eigenfunction. However, this is 
a bad idea since the phase is arbitrary, which can lead 
to strong suppression of the actual eigenfunction with 
respect to numerical errors. It is preferable to remove 
the complex phase factor, using the expression 



:^|3 ^ /(Mec) d^x 
J \p5ec\'^ d^x 



5t 



Sic 



(34) 



It is easy to see that this equation holds for the correct 
eigenfunction, which has constant complex phase. In case 
of numerical errors, the above expression yields an aver- 
aged phase which is insensitive to the unavoidable phase 
errors near the nodes of the eigenfunction, as well as the 
errors at the stellar surface. 

Since analytically the complex phase should be con- 
stant, we can convert it's variation into a measure for 
the quality of the numerical eigenfunction 



d^X 

2j\C\'^d^x ' 



C = pSe. (35) 



It is easy to see that q = for an exact eigenfunction and 
q = 1 in the worst case. Due to the density weight, q is a 
measure for the error of bulk motion, while errors at the 
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surface are ignored. The complex phase of the velocity 
components is removed using the phase computed from 
taking into account the phase shift by 7r/2 of Jw'', 5v^ 
with respect to (5e, Sv'^ . The quality measure is computed 
for each of the velocity components separately. 

Axisymmetric eigenfunctions are extracted from two- 
dimensional simulations, in which case we do not need 
to evaluate Eq. (33). Eigenfunctions of radial modes of 



spherical stars are computed using one-dimensional sim- 
ulations at high resolutions. The eigenfunctions become 
one-dimensional and the integrals ( 34 1 , (351 are replaced 
by integrals over the radial coordinate. 

The eigenfunctions used in this study are only ex- 
tracted once with a good resolution, and then used in all 
simulations, regardless of grid resolution or dimension- 
ality. For this, linear interpolation is used to map the 
numerical eigenfunction to the desired resolution and co- 
ordinate system. This way, we can do convergence tests 
and comparisons between 2D and 3D simulations without 
worrying about differences in the eigenfunctions itself. 

The frequencies are extracted using the time evolution 
of density and velocity at some sample point in the coro- 
tating frame, and of the multipole moments in the inertial 
frame defined in Sec. |II C[ For this, we fit an exponen- 
tially damped sinusoidal, which is usually more accurate 
than using the Fourier transform. The comparison be- 
tween inertial and corotating frame is useful to identify 
m for unknown modes, and serves as a consistency check. 

The frequencies in the inertial and corotating frame 
are related by 



fi l/c - mFn 



(36) 



where fc is the frequency in the corotating frame, fi in 
the inertial frame, and Fn is the rotation rate of the 
star. All frequencies are defined as positive and measured 
with respect to coordinate time. For our setups, Fji and 
fi are also identical to the rotation rate and oscillation 



frequency observed at infinity, see Sec. |IIF| We define m 
such that m > if the wave-patterns of the mode appear 
counterrotating in the corotating frame, and m < for 
modes which are corotating in the corotating frame. Note 
a mode counter-rotating in the corotating frame appears 
corotating in the inertial frame if fc < mFn. 



for models rotating near the Kepler (mass shedding) 
limit. 

For all simulations of a given model, we use one and 
the same spacetime which is computed once with a res- 
olution of at least 100 points per stellar radius. The 
data is mapped onto the computational grid used in our 
simulations using linear interpolation. The shift vector 
is initialized such that we obtain coordinates corotating 
with the (rigidly) rotating star. For three-dimensional 
simulations, we apply the standard transformation from 
cylindrical to Cartesian coordinates. 

To set up spherical stars, we use our own code to solve 
the ordinary differential equations derived by [IB] (TOV- 
equations). For two- or three-dimensional simulations, 
a standard transformation from spherical coordinates to 
cylindrical or Cartesian coordinates is applied. 

We note that the line element found by the two meth- 
ods is not exactly the same for a given nonrotating star 
due to different gauge choices. The line elements for the 
rotating model can be found in [lU . For the simula- 
tion itself the choice of coordinates doesn't matter, since 
it is gauge invariant (up to numerical errors) . It will how- 
ever have an impact on the error of the GW extraction, 
as discussed in Sec. Ill CI 

The time coordinate, which is only fixed up to a global 
factor, is normalized by the initial data codes such that 
a = 1 at infinity. Since the spacetime is stationary, any 
frequency measured with respect to coordinate time at a 
fixed point in the inertial coordinate frame is identical to 
the frequency observed at infinity. 

To excite oscillations, we perturb specific energy e and 
3- velocity w*. Since we are using the Cowling approxima- 
tion during evolution, we neither perturb the metric nor 
reinforce the constraint equations. In case the density 
becomes negative, which can happen near the surface, it 
is simply reset to zero. When perturbing with axisym- 
metric eigenfunctions, we have the freedom to choose the 
phase of the oscillation such that the initial density per- 
turbation is zero, and perturb only the velocity. For non- 
axisymmetric modes, both specific energy and velocity 
are perturbed. 



III. PROBLEM SETUP 



F. Setting up initial data 

Our stellar models are rigidly rotating (and nonrotat- 
ing) stationary configurations of an ideal fluid in general 
relativity, with a polytropic EOS defined by Eq. (12 1. 



The models are characterized uniquely by central density 
Pc rotation rate F^, polytropic exponent F and poly- 
tropic density scale pp. 

To compute the spacetime describing a rigidly rotat- 
ing relativistic star, we use the code described in [44ll45]. 
This code is able to solve the full set of stationary Ein- 
stein and hydrostatic equations with high precision even 



In order to study nonlinear effects, we perturb sta- 
tionary neutron star models with various eigenfunctions, 
which are scaled to amplitudes ranging from the linear 
regime to the strongly nonlinear one, and let the system 
evolve long enough to observe strong damping effects. In 
detail, we study the F, ^fo, ^/2, and ^/_2 modes. 

We use models with three different polytropic EOSs, 
which are summarized in Tab. IH The motivation behind 
our choice is to cover a wide range of stiffness, in order to 
demonstrate it's influence. EOS A was introduced in [47] 
as a rough polytropic approximation to the more realistic 
Pandharipande EOS as tabulated in It is the stiffest 
of our EOSs. The models with EOS B and C (not to be 
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confused with B and C in [JH]) are generic toy models; 
EOS C is very soft, while EOS B is of medium stiffness. 
EOS B is often used as a reference point in numerical 
studies. 



For each EOS, we investigate a nonrotating model 
as well as rigidly rotating models with various rotation 
rates, some close to the Kepler limit. Our models have 
gravitational (ADM) masses 1.4-1.9 Mq, which is in the 
range of observed neutron star masses. We checked that 
the nonrotating models are on the stable branch of the 
mass-radius diagram, and expect the same for the rotat- 
ing ones. The central sound speed is 0.75 c for model 
MAIOO, 0.45 c for MBIOO, and 0.21 c for MClOO. All 
our models are summarized in Tab. HIl 



Model MA 65 is o f particular astrophysical interest. As 
shown in Sec. IV A the counter-rotating mode is most 
probably subject to the CFS instability in full GR. We 
stress that the CFS instability is not active in the Cowl- 
ing approximation. Even in full GR, it's growth timescale 
is on the order of seconds, and therefore irrelevant on the 
timescales of our simulations. However, it could provide 
a mechanism of exciting the high amplitudes investigated 
in this study, provided the instability is not suppressed al- 
ready at much lower amplitudes due to viscosity or mode 
coupling effects. 



TABLE I: The polytropic EOS used in our models, specified 
by polytropic constant V — 1 + 1/n and polytropic density 
scale Pp. The more commonly used polytropic constant K = 



is given in geometric units G - 



Mr.- 



Name 


n 


r 


Pp/e cm"^ 


K 


C 


2 


1.5 


7.000-10" 


2.970 


B 


1 


2 


6.176-10" 


100.0 


A 


0.6849 


2.46 


4.070-10" 


11.65 



TABLE IL Details of the stellar models. M is the gravita- 
tional mass, Fu the rotation rate as observed from infinity, 
Rc the equatorial circumferential radius, Rp and Re the polar 
and equatorial coordinate radius, and pc the central rest mass 
density. 



Name 


EOS 


M/Mq 


Fr/Uz 


Rc/km 


Rp/Re 


Pc/gcm ^ 


MAIOO 


A 


1.615 





9.529 


1 


2.065 


10" 


MA65 


A 


1.910 


1687 


11.61 


0.65 


2.065 


10" 


MBIOO 


B 


1.400 





14.16 


1 


7.905 


lO^-* 


MB85 


B 


1.503 


590.9 


15.38 


0.85 


7.905 


10" 


MB70 


B 


1.627 


792.1 


17.27 


0.70 


7.905 


10" 


MClOO 


C 


1.400 





47.60 


1 


7.618 


10" 


MC95 


C 


1.400 


56.97 


51.17 


0.95 


6.660 


10" 


MC85 


c 


1.400 


83.47 


59.41 


0.85 


5.222 


10" 


MC65 


c 


1.400 


93.11 


80.80 


0.65 


4.007 


10" 



The maximum amplitude of the initial perturbation is 
chosen such that during the evolution, the fluid stays in- 
side a given region, which is usually twice as big as the 
bounding box of the star. Only for 3D simulations of 
model MA65, we had to restrict ourselves to an expan- 
sion factor of 1.6, because otherwise the corners of the 
corotating coordinate system would move with superlu- 
minal speed. 

To excite high amplitude oscillations, we linearly scale 
the eigenf unctions of specific energy and 3-velocity, and 
add them to the background model. For axisymmetric 
simulations, the sign of the perturbation is chosen such 
that the first maximum of the x-velocity along the x- 
axis is positive. For nonaxisymmetric perturbations, the 
sign is irrelevant since changing the sign corresponds to 
a rotation. In the rest of this paper we refer to such a 
high amplitude perturbation based on the eigenfunction 
e.g. of the F mode simply as an _F-mode perturbation. 

Note that this choice is not unique. For example, one 
could scale the momentum density instead of the velocity, 
or change the sign of the perturbation. In the nonlinear 
regime, this will lead to small differences in the results. 
Also, our setups probably differ slightly from what one 
would obtain by letting a mode grow to high amplitudes 
by means of some physical instability or artificial back- 
reaction force. 



IV. NUMERICAL RESULTS 
A. Mode frequencies and eigenfunctions 

As a prerequisite for our studies, we extracted the fre- 
quencies and eigenfunctions of various low-order axisym- 
metric and nonaxisymmetric pressure modes. All fre- 



quencies can be found in Tab. Ill Based on convergence 
tests, including those in |10j, we are confident that the 
numerical error of the frequency is generally below 2 %. 

We also compared our results to those obtained by 
[121 . where a linearized evolution code and Cowling ap- 
proximation was used. The differences are below 1.2 % 
for models MAIOO, MA65, and below 0.6 % for models 
MBIOO, MB85, MB70, which is an excellent agreement. 
For the nonrotating models MClOO and MAIOO, the fre- 
quencies we extracted for the ^/oand modes (which are 
exactly the same analytically) agree better than 0.1 %. 

For the extraction of nonaxisymmetric modes, the 
mode-recycling method is computationally very expen- 
sive, since it requires repeated 3D simulations over many 
oscillation periods. Therefore, we did not repeat the recy- 
cling step as often as for axisymmetric modes, and used a 
resolution of only 50 points per stellar radius, in contrast 
to 100-200 points for 2D simulations. 

As a consequence, when using the nonaxisymmetric 
numeric eigenfunctions as a perturbation, other modes 
are sometimes excited as well, at amplitudes up to a few 
% of the desired mode. We believe that regarding the 
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dynamics of the evolution, this error can be neglected in 
comparison to other numerical errors. When interpreting 
the GW signal on the other hand, it has to be taken into 
account. 

In agreement with [12j . we find that the counter- 
rotating mode of model MA65 becomes corotating in 
the inertial frame. Although we are working in the Cowl- 
ing approximation, this should be the case in full GR as 
well. As shown in [TS] for similar models, dropping the 
Cowling approximation seems to lower the neutral point, 
i.e. the critical rotation rate where the counter-rotating 
/ mode becomes corotating in the inertial frame. There- 
fore, the counter-rotating mode of model MA65 is 
most probably CFS-unstable. 

The cigcnfunction of the mode is shown in Fig- 
ure [l] Although model MA65 is rotating quite rapidly, 
the eigenfunction is only moderately deformed in compar- 
ison to the nonrotating model MAIOO shown in Figure [5] 
The main difference is that, compared to the inner region, 
the oscillation amplitude near the equator is increased. 
This is to be expected since the material is bound less 
strongly with increasing centrifugal force. 

The eigenfunctions we found for stars with moderate 
rotation rates are well described by the slow-rotation 
limit, where the eigenfunctions are given by a spherical 
harmonic times some radial function. For the purpose of 
our study, it is more interesting to look at the rapidly 
rotating models. 

One important observation with regard to GW emis- 
sion is that for rapidly rotating models, the quasiradial 
F mode possesses a considerable quadrupole moment due 
to the oblateness of those stars. As an example, the 
quasiradial F mode of model MA65 is shown in Figure |3] 
We have to stress that the figure shows the eigenfunc- 
tion in the coordinate system set up by the initial data 
code. It would be easy to construct an asymptotically 
Euclidean coordinate system in which even a spherical 
star looks oblate. However, given that MA65 is rotating 
near the Kepler limit, we are confident that the defor- 
mation of the star and the eigenfunction is not mainly a 
coordinate artifact. 

Another effect becomes important very close to the 
Kcplcr limit. As shown in Figure |4j the Vo-niode eigen- 
function of model MC65 is not only strongly deformed 
in comparison to the slow-rotation limit, it also develops 
very pronounced peaks near the equator. Looking at the 
bulk properties of mode, we find that the dynamics of the 
oscillation is still determined by the inner regions of the 
star, but a given amplitude in the interior induces much 
higher amplitudes near the equator. This is not surpris- 
ing since the material there is only marginally bound for 
this model. As a consequence, oscillation modes of MC65 
can store only little energy before nonlinear effects set in, 
as will be shown in Sec. |IVB| In fact, we did not even 
succeed to obtain a clean eigenfunction of the F mode of 
this model, due to strong nonlinear couplings. 

Since the GW luminosity depends strongly on the fre- 
quency, we have to estimate the influence of the Cowling 



TABLE IIL Frequencies fi of the different oscillation modes 
as observed from infinity in the inertial frame, and corre- 
sponding frequencies fc in the corotating frame, q is the 
quality factor of the numerical eigenfunction for 5e defined 
in Eq. ([35|. 



Model 


Mode 


/c/Hz 


/./Hz 


q 




MClOO 


F 


472 


472 


2 


10-'' 


MClOO 


70 


418 


418 


5 


10-'' 


MClOO 


2/. _ 2/. 
72 = 7-2 


418 


418 


4 


10-^ 


MC95 


F 


437 


437 


8 


IQ-'^ 


MC95 


2f 

70 


392 


392 


2 


10-^ 


MC85 


F 


383 


383 


5 


IQ-'^ 


MC85 


2f 

70 


337 


337 


2 


10-^ 


MC85 


2f 
72 


371 


205 


3 


10-^ 


MC65 


2f 

70 


239 


239 


2 


10-^ 


MBIOO 


F 


2686 


2686 


3 


IQ-'^ 


MBIOO 


70 


1883 


1883 


5 


10-" 


Mj30i3 


2-f 

70 


lo93 


lo93 


4 


iU 


MB70 


2f 

70 


1785 


1785 


3 


10-" 


MB 70 


72 


1948 


364 


4 


10-^ 


MAIOO 


F 


4606 


4606 


4 


10-^ 


MAIOO 


2f 

.70 


3038 


3038 


5 


10-" 


MAIOO 


2r _ 2/. 
72 = 7-2 


3037 


3037 


2 


10-^ 


MA65 


F 


3997 


3997 


3 


10-" 


MA65 


2f 

70 


2662 


2662 


2 


10-^ 


MA65 


2f 
72 


2856 


518 


6 


10-^ 


MA65 


2f 

.7-2 


1179 


4553 


2 


10-^ 



approximation. Frequencies in Cowling approximation 
and full GR, computed by means of fully relativistic sim- 
ulations, are given in ^15], for model MBIOO and mod- 
els similar to MAIOO and MA65, amongst others. After 
adding an additional safety margin to account for the 
different models, we assume that we over-estimate the 
frequencies of the F modes by a factor of less than 2.5 
and the ^/o modes by a factor less than 1.5. Further, 
we assume that the frequency of the mode of in the 
corotating frame is over-estimated by a factor < 1.3. For 
model MA65, this implies that the frequency in the iner- 
tial frame is under-estimated by a factor up to 2.3. For 
model MB70, we find that the ^2 mode could become 
CFS-unstable in full GR, with an inertial frame frequency 
in the range 0-90 Hz. 

B. Damping of axisymmetric modes 

For all axisymmetric oscillations of all models, we find 
a common behavior at high amplitudes: if the initial am- 
plitude exceeds a certain threshold, shock waves form 
in the outer layers of the star, which dissipate energy 
during a few oscillation cycles until the oscillation am- 
plitude falls below the threshold again. After this phase, 
we usually observe only the numerical damping. In some 
cases, we also see mode coupling effects, which are how- 
ever small compared to the strongest damping due to 
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FIG. 1: Eigenfunction of the counter-rotating mode of 
model MA65. The left half shows the two-dimensional eigen- 
function 5e{d,z) of specific energy as a contour plot. The 
eigenfunction is zero on the rotation axis. The arrows in the 
right half correspond to the eigenfunction of the velocity in 
the meridional plane. 



FIG. 3: Eigenfunction of the quasiradial F mode of model 
MA65, plotted like in Figure [T] The additional thick solid 
line marks the node. 




FIG. 2: Like Figurejl] but showing the mode of nonrotat- 
ing model MA 100. The specific energy eigenfunction is zero 
only on the rotation axis. 



shocks. 

The final ampUtude can be smaller than the thresh- 
old, in particular for strong initial excitation. The final 
amplitude as a function of the initial one thus has a max- 
imum for some modes. Figure [5] shows a typical evolution 
of the mean velocity in presence of shocks. 

To detect shocks, we monitor the conserved energy Ec 
defined by Eq. ( 30 ) . Figure [6] shows an example of en- 
ergy dissipation in the presence of shocks. To investigate 
the details of shock formation, we produced movies from 
several of our simulations, showing the evolution of the 
density and velocity in the meridional plane. The shock 
formation was clearly visible for simulations with signif- 
icant energy loss. For one-dimensional simulations, we 





FIG. 4; Eigenfunction of the axisymmetric /o mode of model 
MC65. TOP: Eigenfunction Se{d, z) as gray-scale plot. The 
thick solid line marks the node position. The arrows show 
the velocity eigenfunction. BOTTOM: Same eigenfunctions 
scaled by pd, to visualize the bulk properties of the mode. 
The node near the equator is very likely a numerical artifact. 
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tude is independent from the exact value of the evolution 
time. As discussed in Sec. IIIBl the onset of shock for- 



FIG. 5: Evolution of mean velocity v for an _F-mode pertur- 
bation of model MC95 at high amplitudes. 
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FIG. 6: Violation of energy conservation due to shocks, for 
an f-mode perturbation of model MC95. Plotted is the loss 
of conserved energy Ec in units of the total stellar mass versus 
the time in units of the oscillation period To. Only the phase 
of strong decay is shown. 



also verified that shock formation and energy dissipation 
occur simultaneously. The shocks we observed formed 
in the outer layers of the star, but not directly at the 
surface. 

To get an overview on the magnitude of nonlinear 
damping effects, we plot the final amplitude Ap defined 
in Eq. ([28| versus the initial amplitude Aj. As an exam- 
ple, Figure [t] shows the results for the ^/o mode of model 
MC85. 

In general, we choose different evolution times for dif- 
ferent sequences, such that the strong damping phase at 
high amplitudes (compare Figure [5| is included, but not 
much longer. We choose short evolution times because 
we are interested mainly in strong damping effects. In 
any case, we can only accurately measure damping ef- 
fects acting on timescales significantly shorter than the 
timescale of numerical damping. 

For each sequence, we extract the amplitude Ao &t 
which nonlinear damping effects become stronger than 



mation is predicted accurately despite the use of a cold 
(one-parametric) EOS. Therefore, the onset of nonlinear 
damping is not affected by this approximation neither. 

The values we found for the Fand ^/o modes of the dif- 
ferent models cover a wide range Ad ~ 0.003 . . . 0.13 c. 
We will discuss the dependence on the model parame- 
ters in sections IV C and IV D| Naturally, the values of 
Ap, provide upper limits for the applicability of linear 
perturbation theory. 

The final amplitude contains only information on how 
much the system is damped, but not how fast. To quan- 
tify the damping speed, we use the timescale r of the 
initial decay defined in Eq. (29 1. Figure [s] shows r versus 
the initial amplitude Aj for the case of the axisymmetric 
^/o mode of model MA65. We stress that the values of 
T are affected by the use of a cold EOS, as discussed in 



the numerical damping, given in Tab. IV This ampli- 



Sec. |IIB[ and should be regarded as an educated guess. 

Not surprising, the damping becomes faster with in- 
creasing amplitude. By definition, r is a measure for time 
averaged decay. More detailed investigation reveals that 
for the highest amplitudes, the larger part of the energy 
is dissipated on timescales smaller than the oscillation 
period during the formation of a strong, but short-lived 
shock. 

The decay timescale can be used to estimate the re- 
quired growth timescale of a hypothetical instability sat- 
urating at a given amplitude. We are not aware of any ax- 
isymmetric instability acting on timescales shorter than 
the numerical damping. For any process with longer 
growth times, the value Ad provides an upper limit for 
the saturation amplitudes. 

To estimate the numerical errors, we computed sev- 
eral sequences using 3 different resolutions. Note it is 
important to choose the sign of the perturbation con- 
sistently, since we observed differences in the nonlinear 
regime comparable to the numerical error at low resolu- 
tion. 

By far the greatest errors are found for the stiff EOS A. 
Figure [S] shows the convergence of the nonlinear damping 
timescale of the •^/g mode of model MA65 as a represen- 
tative example. To explain the behavior at low ampli- 
tudes, we note that when approaching the linear regime, 
the timescale of physical damping goes to infinity. In 
that case, we only see the numerical damping, which be- 
comes weaker with increasing resolution. The stronger 
the damping, the more accurate is the value we find for 
r. 

For the soft EOS C, the numerical damping is ex- 
tremely low. For a resolution of 100, we find timescales 
tm > 800 ms in the linear regime. Also the nonlinear 
results match extremely well, as shown in Figure [7] The 
accuracy of models with EOS B is in between the other 
two cases. Since the numerical damping is mainly given 
by surface effects, as shown in [10], and the steepness of 
the density gradient at the surface increases with stiff- 
ness, this behavior is not surprising. 
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FIG. 7: Nonlinear damping of the ^/o mode of model MC85. 
Plotted is the final amplitude Af after an evolution time of 
20 ms over the initial amplitude Ai. We show results of ax- 
isymmetric simulations with resolutions of A'^ = 50, 100, and 
200 points per stellar radius, as well as a three-dimensional 
simulation with resolution A'' = 50. 



TABLE IV: Onset of nonlinear damping. An is the initial 
amplitude Ai at which nonlinear damping is strong enough 
to be distinguishable from numerical damping, read off from 
plots of initial versus final amplitude. 



Model 


Mode 


Ab/(10-3 c) 


MClOO 


F 


11 ± 1 




% 


4± 1 


MC95 


F 


9± 2 




2-f 

70 


5±2 


MOoo 


r 


7 ± 2 




2-f 

70 


4± 1 


MC65 


2f 

70 


< O.l 


MBIOO 


F 


53 ±7 




2-f 

70 


20 ±5 


MAIOO 


F 


100 ± 25 




70 


45± 15 


MA65 


F 


40± 15 




2-f 

70 


8±3 
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FIG. 8: Timescale of initial decay t due to shock damping 
versus the initial amplitude Ai, for an ^/o-mode perturba- 
tion of model MA65. Shown are results from axisymmetric 
simulations with three different resolutions as well as one 3D- 
simulation of the same problem. The resolution A'^ is given in 
points per stellar equatorial radius. 



To get a rough error estimate for siraulations which we 
only ran at one resolution, we assume that the numerical 
damping observed in the linear regime for single oscil- 
lations is also operational at high amplitudes, acting on 
the same timescale tjj. We then scale our time series by 
exp(t/TD) to obtain corrected values for final amplitude 
Ap and damping time r. 

To estimate the accuracy of three-dimensional simula- 
tions we also performed 3D simulation of axisymmetric 
modes at high amplitudes. We generally found a good 
agreement with the axisymmetric results for the same 
setup; see Figure [7] and Figure [8] 



C. Influence of the EOS on damping 

To study the influence of the EOS, we use the nonrotat- 
ing models MClOO, MBIOO, MAlOO, which have similar 
masses, but EOSs of different stiffness. 

The amount of nonlinear damping for the different 
models is shown in Figure [9j The amplitudes Au at 
which nonlinear damping effects become visible is given 
in Tab. Jvl 

As one can see, the nonlinearity sets in one order of 
magnitude earlier for the EOS C, which is the softest 
one, than for EOS A, the stiffest one. Also the max- 
imum final amplitude is the largest for the stiff EOS. 
Note that MAIOO is 13 % heavier than MClOO, which 
might account for part of the differences. Nevertheless, 
the masses of models MClOO and MBIOO are exactly the 
same, and still the damping differs strongly. We conclude 
that with increasing stiffness of the EOS a neutron star 
(of a given mass) can store more energy in axisymmetric 
oscillations before strong damping sets in. 

We believe that the onset of shock formation is deter- 
mined mostly by the stiffness of the EOS in the density 
range of the outer layers of the star, since this is the 
region where the shocks form. However, this requires 
further investigation. 

In Figure |10[ we plot the damping timescale for 
the different nonrotating models. At a given damping 
timescale, the possible oscillation amplitudes differ by 
roughly one order of magnitude. Further, the damping 
timescale at the highest possible amplitudes (before ma- 
terial starts leaving the computational domain) decreases 
with increasing stiffness, i.e. the damping is faster. 
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FIG. 9: Onset of nonlinearity for the _Fand /o modes of the 
nonrotating models MBIOO, MClOO, and MAIOO. Plotted is 
the final amplitude Ap versus the initial one Ai. The solid 
gray lines are values corrected for the numerical damping. 
The F mode results are extracted from one-dimensional sim- 
ulations with a resolution of 200 points per stellar radius, the 
^/o modes are computed with two-dimensional simulations at 
resolution of 100^. 
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FIG. 10: Initial damping timescale t versus initial amplitude 
j4/, for the same simulations shown in Figure [9] The solid 
gray lines are values corrected for the numerical damping. 
We only show amplitudes for which the difference is less than 
50 %. For slower damping, the error of r blows up quickly. 



D. Influence of rotation on damping 

To investigate the influence of rotation on axisymmet- 
ric modes, we first compare the models MClOO, MC95, 
and MC85, which have same mass and EOS, but differ- 
ent rotation rates up to 90 % of the Kepler limit. The 
amount of damping for the Fand ^/o modes is shown in 
Figure [TT] and the damping timescales in Figure [T2j For 
the stiff EOS A, we compared nonrotating model MAIOO 
with the rapidly rotating model MA65. Note however 
that MA65 is also 18 % heavier. The results are shown 
in Figure [T3| 

The differences between models MClOO and MC95 are 
very small. For MC85 on the other hand, the nonlinear 
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FIG. 11: Influence of rotation on nonlinear damping of ax- 
isymmetric modes. Shown are final versus initial amplitude 
for a perturbation with the F mode (TOP) and the '^/o mode 
(BOTTOM) of models MClOO, MC95, MC85. The evolution 
time was 20 ms for the ^/o mode and 15 ms for the F mode. 



damping is significantly faster and stronger in compar- 
ison. The onset of nonlinearity occurs at slightly lower 
amplitudes for the F mode, while for the ^/o modes we 
found no significant difference. Since the fluid is bound 
less strongly with increasing centrifugal force, an increase 
of nonlinearity at given amplitude is to be expected. 

An extreme case is given by model MC65, which is 
rotating almost at breakup velocity. For this case non- 
linear effects start much earlier than for MC85. For 
the ^/o mode, nonlinear effects are definitely present for 
Ai > 10"''. At an amplitude as low as Aj = 1.5 • 10~^ a 
fraction of 2 • 10^'"' of the total mass leaves the computa- 
tional domain. We obviously observe what is called mass 
shedding induced damping, i.e. even small oscillations 
liberate mass near the equator. The same mechanism 
has already been observed in |H1 1^ for models with EOS 
B. Since stars rotating so close to the Kepler limit are 
unlikely to occur in nature, we did not investigate this 
case further. 



E. Damping of nonaxisymmetric modes 

To study the evolution of nonaxisymmetric modes, we 
performed three-dimensional simulations of the counter- 
rotating mode of the rotating models MA65, MAIOO, 
MB70, and MC85. For model MA65, we also studied the 
corotating ^/_2 mode. 

In contrast to the axisymmetric case, we did not ob- 
serve any steep decrease of the conserved energy that 
would point to formation of strong shocks. We did how- 
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FIG. 12: Influence of rotation on damping timescale r of 
F modes (TOP) and ^/o modes (BOTTOM). The values are 
already corrected for the numerical damping. Only points are 
shown for which the correction is less than 50 %. 
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FIG. 13: Nonlinear damping of F mode (TOP) and mode 
(BOTTOM) for nonrotating model MAIOO and rapidly rotat- 
ing model MA65. Evolution time was 4 ms. 



ble strength at the resolutions we could afford. 

To estimate the error of the damping, we assume that 
the timescale of the numerical damping in the nonlinear 
regime is equal to the one observed in the linear regime. 
For the axisymmetric case described in Sec. IV B this 
estimate agreed well with the convergence test results. 
For axisymmetric setups, we also proved that the code is 
as accurate in 3D as in 2D. As a further check, we com- 
pared resolutions of 50 and 75 points per stellar radius. 
The results are shown in Figure [14] to [16] 

To identify the damping mechanism, we produced an- 
imations showing the density and velocity evolution in 
the equatorial plane as well as the meridional planes. 
For high amplitudes, we found that the wave-patterns 
traveling around the star become more and more non- 
sinusoidal towards the surface, where we observe effects 
similar to wave breaking in ocean waves. A snapshot is 
shown in Figure [18] and [19] We do not find any formation 
of shocks in the interior of the star, only at the surface 
in regions with densities around 10^'' . . . lO^'^ times the 
central one. Note this also justifies the approximation 
of adiabatic evolution implied by the use of a polytropic 
EOS. For low amplitudes, the wave-patterns looked ex- 
actly like the eigenfunctions we used to excite the os- 
cillation. One important dissipative mechanism at high 
amplitudes thus seems to be wave breaking at the sur- 
face. As will be shown in the next section, the ^2 modes 
can also participate in mode coupling. 

Since wave breaking is very difficult to resolve numeri- 
cally, our results are probably not as accurate as the ax- 
isymmetric ones. More important, the surface of realistic 
neutron stars does not consist of a boundary between a 
cold fluid and vacuum, like in our models; there is either 
a solid crust or a hot envelope, and also a magnetic field. 
To compute the damping of the nonaxisymmetric modes 
correctly, it will be necessary to add a more detailed de- 
scription of the surface to the numerical models. 

Taking into account the strong deformation of the star, 
one also has to question the accuracy of the Cowling ap- 
proximation. On the other hand, the gravitational field 
is mainly determined by the denser regions of the star, 
for which we observe much smaller deformations. We 
therefore assume that the error on the bulk motion of 
the star is comparable to the linear regime. This should 
be verified by means of a fully relativistic study. 



F. Mode coupling 



ever observe a continuous decrease of the mean velocity 
and the I = m = 2 multipole moment. The dissipation 
of energy was considerably slower than the one observed 
in case of shock formation in axisymmetric oscillations, 
even for high perturbation amplitudes. 

Figures [T4] [T7] show the damping timescale. For mod- 
els MA65 and MB70, it proved difficult to disentangle 
physical and numerical damping, which are of compara- 



For the % modes of models MA65 and MB70, we 
found that beside wave breaking, mode coupling has a 
significant influence on the damping. 

In the spectrogram shown in Figure |20[ one can eas- 
ily spot mode coupling between the mode and some 
unidentified mode, called mode X in the following. For 
low amplitudes, only the mode is excited. With 
increasing amplitude, mode X is growing faster. At 
the highest amplitude, it becomes dominant in the v'^- 
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FIG. 14: Timescale r of initial damping for the /2 mode of 
model MC85. Shown are numerical results obtained at resolu- 
tions 50 and 75 points per stellar equatorial radius. The error 
bars are our estimate for the influence of numerical damping. 
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FIG. 16: Like Figure 14 but for the 72 mode of model MA65. 
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FIG. 15 

MB70. Error bars leaving the plot are of infinite size, i.e 



but for the 



mode of model 
the 



numerical results are compatible with no physical damping. 



spectrum on a timescale as low as 4 ms. 

The Fomier spectrum for amplitude Aj = 0.04 c is 
Mode X is located at 1458 ± 50 Hz, 



shown in Figure 21 



the F mode at 3997±50 Hz, the 7o mode at 2590±50 Hz, 
and the mode at 2912 ± 50 Hz. 

From this, we find two interesting resonances for the 
unknown mode. First, mode X is located at half the 
frequency of the mode, measured in the corotating 
frame. Second, the frequency also equals (up to the res- 
olution of the Fourier spectrum) the difference between 
the _Fand the ^/o mode. The latter could be interpreted 
as a hint that the unidentified peak does not belong to 
a proper oscillation mode, but is a nonlinear harmonic 
instead, i.e. a typical feature of Fourier transforms in 
case of nonharmonic waveforms and nonlinear superpo- 
sitions. However, the F mode is not present in the spec- 
trum of w"^, and the ^/o mode is neither present in the 
spectrum of v'^ nor of w'', while the unknown peak is 
present in the spectra of v*^ and v"^ . This supports the 
interpretation that the peak belongs to an actual oscil- 
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FIG. 17: Nonlinear damping of nonaxisymmetric modes for 
nonrotating model MAIOO and rapidly rotating model MA65. 



lation mode, because the amplitude of a nonlinear har- 
monic in the Fourier spectrum of some quantity scales 
with the product of the amplitudes (in the same spec- 
trum) of the modes producing the harmonic. However, 
it is still possible that the unidentified peak in the spec- 
trum of the multipole g20 is a nonlinear harmonic of the 
Fand ^/q modes, and thus unrelated to mode X. 

Since the frequency of the unidentified mode is in 
the inertial mode range, and well below the frequencies 
of the lowest order axisymmetric pressure modes, it is 
most probably an inertial mode. In order to identify 
the modes appearing for models MB70 and MA65, we 
first attempted mode recycling with axisymmetric sim- 
ulations. However, the method did not converge. This 
does not necessarily mean that the mode is nonaxisym- 
metric. As described in jl 1] . extracting inertial modes 
using mode recycling is more difficult than for pressure 
modes. Nevertheless, several low-order inertial modes of 
model MB70 were identified in [TT] ; none of them matches 
the frequency of the unknown mode. Since extracting in- 
ertial modes typically requires more recycling steps than 
pressure modes (due to the dense spectrum), and because 
the frequencies are lower, extracting nonaxisymmetric in- 
ertial modes is computationally expensive. Therefore, we 
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FIG. 18: Snapshot of the evolution of density and velocity in 
the equatorial plane, for an ^/2-mode perturbation of model 
MB70 with amplitude Aj = 0.04. Shown are the lines of 
constant density p, spaced at regular intervals of y/p. The 
outermost line corresponds to a density of 10^* times the 
central density. The arrows correspond to the velocity in the 
corotating frame. The wave-patterns are rotating clockwise, 
leaving behind low density regions which fall back. Note those 
regions are not well resolved numerically. 



only performed the first step, using an ^/2-niode pertur- 
bation of model MA65 at amplitude Aj = 0.04, for which 
we know that the unknown mode is strongly exited. The 
resulting three-dimensional numerical eigenfunction had 
predominantly a \m\ — 2 0-dependency. We tried ex- 
tracting a two-dimensional eigenfunction assuming either 
TO = +2 or TO = —2, but obtained a big phase error 
{q ~ 0.5), which typically points to a superposition of dif- 
ferent modes of almost the same frequency. In our case, 
it seems that the unidentified peak actually consists of 
two high order inertial modes with m = ±2, although a 
clear identification will require additional mode recycling 
steps. 

We note that the mode coupling is significantly 
stronger for resolution 75 than for resolution 50, which 
means the process is not well resolved numerically. In 
the aforementioned low-quality numerical eigenfunctions, 
we found structures smaller than one fifth of the stellar 
radius. Therefore, the excited inertial modes probably 
suffer from strong numerical damping, which suppresses 
the coupling at resolution TV = 50. Further, the time 
evolution in presence of mode coupling effects typically 
depends strongly on the initial amplitudes of both modes, 
even if one amplitude is small. In our case, the unknown 
mode is probably seeded only by numerical errors; in real 
stars, the situation might be different. It is hence possi- 
ble that the mode decays even faster in reality than 
in our simulations. 

Beside the inertial modes, the mode seems to cou- 




FIG. 19: TOP: Like Figure [l8| but for the mode of model 
MA65. BOTTOM: Same data, but only densities above 10"^ 
times the central one are shown, and the arrows correspond 
to the momentum density instead of velocity. 



pie to the V-2 mode as well, albeit less dramatic. This 



coupling will be discussed in Sec. |IVG[ since it interferes 
with the estimation of the GW strain. It is still unclear 
what the exact nature of the observed couplings is, and 
whether the mode coupling picture is adequate at such 
high amplitudes at all. 

To measure the decay in presence of mode coupling, v is 
not useful because it is quite insensitive to transfer of en- 
ergy between modes. However, we can use the dominant 
multipole moment of the mode, which for the modes 
is the I — m = 2 mass multipole. In the Fourier spectra 
of (722, the dominant contribution is due to the ^/2 mode. 
By fitting a function q'^ exp{—t/Tm + ii^t) to the initial 
evolution of the multipole moment ^22, we obtain another 
damping timescale r^. 

While T is a good measure for the overall energy dis- 
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FIG. 20: Mode coupling at different amplitudes of the 
mode of model MA65 to an unknown mode at half the 
frequency. Shown is time evolution of the Fourier spectrum 
of the 0- velocity in the corotating frame, at r = |i? on the 
space diagonal. The spectra were computed using a running 
window of width 4 ms. The perturbation amplitudes Aj are 
0.01 (upper left), 0.02 (upper right), 0.04 (lower left), and 
0.08 (lower right). The frequency resolution is 250 Hz. 
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FIG. 22: Decay timescales r and Tm of v and q22 versus initial 
amplitude Aj, for models MB70 and MA65 perturbed with 
the mode. The results are obtained at resolution 75. The 
error bars estimate the influence of the numerical damping. 



TABLE V: Upper limits on the onset of nonlinearity, for non- 
axisymmetric oscillation modes. An is the initial amplitude 
Ai at which nonlinear damping (including mode coupling ef- 
fects) can be safely distinguished from the numerical damp- 
ing, r is an upper limit for the decay timescale of the mode 
at amplitude Ad- 



Model 



Mode 



An/{W-'' c) 



r/ms 



MC85 
MB 70 
MAIOO 
MA65 
MA65 



^2 

J2 

7-2 



10 

35 

130 

18 

30 



90 
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20 
20 

40 



star, this is only possible if 



FIG. 21: Fourier spectra of u*^, , and 520, in arbitrary 
units, for a ^/2-mode oscillation of model MA65, at amplitude 
Aj = 0.04. All frequencies are given in the corotating frame. 



sipation, r„ is a measure for the decay of the mode. 
Figure 22 shows a comparison between the two. As one 



(37) 



can see, mode coupling is even more important for the 
decay of the ^2 mode than dissipative effects, e.g. wave 
breaking. 

At least this is the case for models MA65 and MB70; 
for the mode of model MAlOO and the mode of 
model MA65, we do not observe mode coupling, while 
for model MC85 we observe weak coupling to a mode at 
one quarter of the (corotating) frequency of the ^/2 mode. 
What is the difference? In model MB70 and MA65, the 
coupled mode is an inertial mode which oscillates at ^/c, 
with /c being the frequency of the ^/2 mode in the coro- 
tating frame. Since the inertial mode range is bounded 
by 2Fr (see [H]), where is the rotation rate of the 



Model MC85 (and MAlOO anyway) rotates too slowly to 
satisfy this condition. The 1:4 resonance to an inertial 
mode on the other hand is possible for model MC85. 

Note the above condition automatically holds for CFS- 
unstable ^2 modes, for which fc < 2Ffl. It seems plau- 
sible that the coupling is also present at lower ampli- 
tudes, but evolves at longer timescales. To compute this 
timescale and thus a more stringent limit on the satura- 
tion amplitude of CFS-unstable ^/2 modes, perturbative 
studies are required. 

The upper limits for the onset of nonlinearity and the 
corresponding damping timescales for the nonaxisym- 
metric modes are summarized in Tab. [V] The mode cou- 
pling effects for models MA65 and MB70 are taken into 
account. Since the CFS instability of the ^/2 mode of 
model MA65 has a growth time in the order of seconds, 
our result yields an upper limit for the saturation ampli- 
tude. 
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G. Gravitational luminosity 

To estimate the GW production in our simulations, we 
use the multipole formahsm described in Sec. |II C| For 
each simulation, we plot the time evolution of the mag- 
nitude of the strain coefhcients Af^ and Af^ , as shown 
in Figure [23j 

In general, the only significant contribution is due to 
the I = 2 mass multipoles 522 or (720- Unless noted other- 
wise, one of the strain amplitudes A2Q or A22 is strongly 
dominant, and will be denoted by A in the following. For 
the axisymmetric case, terms with m 7^ vanish. 

For oscillations in the linear regime, the strain is pro- 
portional to the amplitude. We introduce a normalized 
mode strain A and luminosity C by 



A 



10^ 



-A, 



L 



10- 



C. 



(38) 



From the simulations with the lowest amplitudes, which 
are in the linear regime, we extracted A and C for each 



mode. The results are given in Tab. VI 

The strain produced by the ^/2 mode of model MA65 
is the lowest one of all modes. This is mainly due to 
the low frequency in the inertial frame. The corotating 
^/_2 mode produces a strain larger by a factor of 143 at 
the same amplitude; the frequency dependency ac- 
counts for a factor of 77. 

Since model MA65 is already close to the Kepler limit, 
the strain cannot be increased significantly by uniformly 
spinning up the star. It might however be possible to con- 
struct differentially rotating models for which the CFS- 
unstable ^/2 modes have higher frequencies. 

The strain of the mode of model MC85 is surpris- 
ingly high given that it's frequency is even lower than for 
model MA65. On the other hand, the radius of model 
MC85 is 5.1 times the one of MA65, which enhances 
the I — 2 multipole moments by a factor 25. We thus 
consider the high strain amplitude of model MC85 as a 
consequence of the unrealistically large radius. 

Another noteworthy result is that the strain ampli- 
tudes of F modes of rapidly rotating stars are c ompa rable 
to those of the ^/o modes. As discussed in Sec. II C can- 
cellations in the integrands of the multipoles prevented 
us from computing a meaningful error estimate for the 
strain produced by the F modes. We nevertheless believe 
that at least for model MA65, the high strain amplitude 
of the F mode is not just an artifact due to the use of the 
quadrupole formula, but a genuine effect of the rapid ro- 
tation rate. However, this can only be verified by means 
of a fully relativistic study, using direct wave extraction 
methods. The question is relevant for GW astronomy 
because strong _F-mode oscillations are likely to occur in 
a number of astrophysical scenarios. 

An overview over the maximum GW strain observed in 
the nonlinear simulations of axisymmetric modes can be 
found in Figure [24} As one can see, the maximum strain 
roughly scales with the amplitude; by simply extrapo- 
lating from the linear regime, one obtains an estimate 
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FIG. 23: Evolution of GW strain amplitudes for model MA65 
perturbed with the corotating mode at amplitude Ai = 
0.01 c. Shown are only the 5 biggest contributions. 



TABLE VI: Gravitational luminosity and strain produced 
by oscillation modes in the linear regime, normalized to an 
amplitude v — 10"'^. The values are computed from a fit to 
the multipole moments. The err or bo unds due to the use of 

the multipole formalism (see Sec. II C I are Aq~^ . . . Aq for the 

-2 r„2 



strain, and Cq . . . Cq for the luminosity. For the F modes, 
the error is unknown; the values are given only for complete- 



ness. 


Model 


Mode 


A/10 Mpc 
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MA65 
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2.96 


lQ-25 


1.06- 
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7 


MA65 
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1.24 


lQ-24 


8.28- 
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13 


MA65 
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4.09 
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1.37- 
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13 


MA65 




5.84 
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MAIOO 
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3.59- 


10*2 
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MAIOO 


2f 
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4.96 
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6.85- 


10*2 
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MBIOO 
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7.79 
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1.63- 


10*2 


1.5 


MB70 


2f 
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1.10 
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MClOO 


2f 

Jo 


5.28 
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MClOO 


2f 
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3.74 


10-24 


7.39- 


1040 


1.2 


MC95 


F 


6.11 


10-25 


5.38- 


1038 


? 


MC95 


2f 
JO 


4.08 


10-24 


1.93- 


1040 


1.5 


MC85 


F 


1.79 


10-24 


3.54- 


1039 


? 


MC85 


2f 

Jo 


3.41 


10-24 


9.95- 


1039 


1.4 


MC85 


2f 
J2 


8.65 


10-25 


9.44- 


1038 


1.4 



which is correct up to a factor of 2 even for the highest 
amplitudes. 

For each simulation, we computed the Fourier spectra 
of the strain to determine the frequency of the dominant 
contribution. Usually, the only significant contribution is 
due to the mode used for perturbation. For the mode 
of MA65 on the other hand, we find that the GW strain 
is mainly due to the corotating ^/_2 mode, as well as the 
axisymmetric Fand ^/o modes. Those modes have a lower 
amplitude than the ^/2 mode. But since the mode of 
this model is a relatively inefficient GW emitter, due to 
the low oscillation frequency, the other modes dominate 
the GW strain. 
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FIG. 24: Maximum gravitational wave amplitude A = \A2o\ 
versus initial amplitude Aj, for axisymmetric perturbations. 
The contribution of other multipole moments to the strain are 
negligible in all cases. The dotted lines are an extrapolation 
from the lowest amplitude, assuming linearity. 



FIG. 25: Spectrogram of the gravitational wave amplitude, 
for the ■^/2-mode oscillation of model MA65, at perturbation 
amplitudes Ai = 0.01 (TOP) and Ai = 0.04 (BOTTOM). 
Shown is the Fourier spectrum of A20 (LEFT) and 5R(Ai?) 
(RIGHT) inside a running time window of width 4 ms, versus 
the time at the center of the window. 



To understand the strain of the mode, we thus need 
to understand why the other modes are present. Fig- 
ure [25] shows a spectrogram of the corresponding strain 
amplitudes. By further comparing the absolute ampli- 
tudes of the multipole moments ^22 and (720, we find the 
following picture: 

The mode is excited due to a small contamination 
of the numerical eigenfunction; it's initial amplitude is 
around 4 % of the ^/2 mode in the linear regime. For 
high amplitudes however, it grows up to 20 % during the 
evolution. The reason is unclear, it might be another 
mode coupling effect. 

The F- and ^/o modes are present from the beginning 
as well, but in the linear regime their amplitude becomes 
negligible compared to the amplitude of the ^/2 mode. 
For a high initial amplitude Aj = 0.04 c, the amplitude of 
(720 is still less than 20 % of (722 • We explain the presence 
of the .Fand ^/o modes by the ad hoc way in which we 
scale the eigenfunctions to high amplitudes; it would be 
surprising if only one mode was excited. The relative 
amplitudes could be totally different for the case of a 
^/2 mode slowly grown to high amplitudes by the CFS 
mechanism. 

To quantify the strain produced by the ^/2 mode alone, 
we fit a function qg exp(— i/r^ -\-iijjt) to the (722 multipole 
and then compute the strain from the multipole moment 
given by the fit. This way, we effectively filter out contri- 
butions of frequencies other than the one of the ^/2 mode. 
Note all the strains in the linear regime, given in Tab. |VI[ 
have been computed in the same way (for the axisym- 
metric modes, we use a damped sinusoidal since (720 is 
real- valued) . 

The results for different amplitudes are shown in Fig- 
Only 3-24% of the total strain is due to the main 



ure 
2 



26 



■/2 mode. We stress that the total strain found for our se- 
tups is probably higher than the strain of a ^/2-oscillation 
grown to high amplitude by the CFS-mechanism. To ac- 
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FIG. 26: Initial GW amplitude A versus initial oscillation 
amplitude Aj, for model MA65 perturbed with the counter- 
rotating mode. Shown is the total strain as well as the 



strain due to the /2 mode alone (see main text) 
lines are extrapolations from the linear regime. 



The dotted 



curately compute the latter, one needs to model carefully 
the growth of secondary coupled modes, even if they are 
dynamically unimportant. 

For the V2 mode of model MB70, we qualitatively find 
the same picture as for model MA65. For the V2 mode 
of models MC85 and MAlOO, as weU as the V_2 mode 
of MA65, contribution of other modes to the GW strain 
are small. The strain at high amplitudes is shown in Fig- 
ure |27[ Again, the extrapolation from the linear regime 
yields a rough estimate, correct up to a factor of 3. 



H. Gravitational wave detectability 

In the following, we estimate the prospects of detecting 
GW in the case of either a one-time strong excitation or 
a continuous excitation at lower amplitudes. 
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FIG. 27: Maximum gravitational wave amplitude A - 
versus initial amplitude Ai, for m = 2 perturbations. For 
the ■^/2 modes of models MB70 and MA65, only the contri- 
bution of the ■^/2 mode itself is plotted (see main text). The 
dotted lines are an extrapolation from the lowest amplitude, 
assuming linearity. 



First we investigate the continuous oscillation of one 
mode at saturation amplitude of some instability. Our 
primary interest here is the CFS instability of the 
mode of model MA65. This instability has a growth 
time of seconds, which is longer than the timescale of the 
numerical damping. Therefore, we can only give upper 
limits on the saturation amplitude, namely the ampli- 
tude at which we definitely see nonlinear damping. The 
actual saturation amplitude may be much smaller. 

Figure 28 shows the gravitational strain at the ob- 



served onset of nonlinear damping for several oscillation 
modes. We plot the hee component given by Eq. (14), 
assuming an optimal viewing angle. The general expres- 
sion for the spin tensor harmonics can be found in |49| 
(note the standard notation of magnetic and electric type 
radiation is reversed there). 

The detectability of sinusoidal GW signals depends not 
only on the strain amplitude, but also on the integration 
time. Effectively, the detectable amplitude scales with 
•\/]V, where N is the number of available wave cycles. 
For a continuous signal, this only depends on the detec- 
tor and the computational resources for data analysis. If 
we assume at least a few hundred cycles, the mode at 
the onset of strong damping is detectable with advanced 
LIGO at a distance of 10 Mpc. Under the assumption 
that the saturation amplitude of the CFS instability of 
the mode is only limited by the strong damping ef- 
fects investigated in this work, it would be possible to 
detect sources as far away as the Virgo galaxy cluster. To 
compute the actual saturation amplitude however, per- 
turbative studies at lower amplitudes are needed. 

The numbers above are valid for model MA65, which 
is an extreme case. Models rotating slightly slower will 
produce considerably lower strain. The reason is that the 
strain scales with the square of the mode frequency in the 
inertial frame, which for the modes crosses zero when 



the rotation rate falls below the neutral point. Keeping 
mass, EOS, and mode amplitude v fixed, we can make an 
expansion of the strain amplitude A in terms of rotation 
rate near the neutral point: 



A = Ao 



(39) 



where Aq, are the strain amplitude and rotation rate 
of MA65, and F^ is the rotation rate at the neutral 
point. Further, it is safe to assume that the amplitude 
marking the onset of nonlinearity behaves smooth near 
the neutral point, since the rotation rate influences the 
dissipation effects mainly via the centrifugal force, and 
the inertial modes participating in mode coupling depend 
on the Coriolis force; both forces depend on the absolute 
rotation rate. 

For frequencies in the range 0.1-2 times the ^/2-mode 
frequency of model MA65, the sensitivity curve of ad- 
vanced LIGO is quite flat. In this range, the upper limit 
we find for the distance at which we can observe the 
mode roughly scales with (Fr — F^)"^. 

Since model MA65 is close to the Kepler limit, the 
strain cannot be increased significantly by increasing the 
rotation rate; on the contrary, we expect that nonlinear 
effects set in much earlier when approaching the Kepler 
limit further, as observed for model MC65. Further, we 
doubt that one can construct models where the frequency 
of the CFS-unstable ^/2 modes is more than three times 
higher without resorting to extremely unrealistic EOSs. 
On the other hand, this statement only applies to uni- 
formly rotating stars; it might be fruitful to investigate 
differentially rotating models as well. 

Another restriction for detection of continuous signals 
is currently given by the limited computing resources for 
data analysis, which prevents to search the whole pa- 
rameter space for signals. In particular, searches [50} [51] 
for continuous signals in early data from the 5th sci- 
ence run of LIGO assumed a maximum time derivative 
of the signal frequency. Those searches where targeted 
to signals from spinning nonaxisymmetric neutron stars, 
based on estimates on the maximum sustainable elliptic- 
ity. The luminosity produced by such sources is orders 
of magnitude below the upper limit we computed for the 
CFS-unstable ^/2 mode of model MA65. Accordingly, 
our models spin down much faster. Assuming that the 
mode frequency fc in the corotating frame stays con- 
stant, the decrease of the signal frequency, i.e. the mode 
frequency in the inertial frame, is fi = mFfi, provided 
that fc < mFfi. The spin-down rate can be computed 
from the expressions for the angular momentum lost by 
gravitational waves given in [42^ and the moment of iner- 
tia computed by the initial data code. For the ^/2 mode 
of model MA65 at the upper limit, wc compute a sig- 
nal frequency decrease of / = —5-10^^ Hz/s, which 
is orders of magnitude outside the ranges searched in 
[5CTI 1^ . For models closer to the neutral point, lumi- 
nosity and spin-down rate would be smaller. As a rough 
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FIG. 28: Gravitational wave strain hee and mode frequency, 
for ^/2and modes at amplitudes where we are confident to 
observe nonlinear damping. The strain is given for a source 
at distance 10 Mpc, assuming optimal viewing angle. The 
lines are estimates for the error due to the Cowling approxi- 
mation, assuming that only the frequency changes. The error 
bars correspond to the error due to the use of the quadrupole 
formula. Other errors are negligible. For comparison, we plot 
the sensitivity curves of advanced LIGO and the proposed 
Einstein Telescope (ET). 



estimate, we assume that only the frequencies change, 
while maximum amplitude and oscillation pattern stay 
constant. The spin-down rate then scales with ff. It is 
still above the maximum for the whole frequency range 
50 . . . 1100 Hz investigated in [50l[5T]- Those searches are 
clearly only sensitive to neutron star oscillations with am- 
plitudes well below the onset of strong damping effects 
investigated here. 

For modes other than the mode of model MA65, 
we are not aware of any instability. Nevertheless, we pro- 
vide limits for a hypothetical instability. Figure [28] and 
29 shows upper limits for the strain amplitudes assuming 
that the instability acts on timescales slower than the nu- 
merical damping. Figure [30| gives the strain at saturation 
amplitude for the optimistic assumption that the hypo- 
thetical instability has growth time t/ = 20 Tp, where Tp 
is the period of the F mode of the corresponding stellar 
model. 

Besides slowly growing instabilities, oscillations can be 
excited by sudden violent events, e.g. glitches, magnetar 
giant flares, phase-change induced collapse ([55]), and 
mergers forming a metastable rapidly rotating star. It is 
beyond the scope of this article to estimate the ampli- 
tude of oscillations excited by such events. We investi- 
gate the detectability for a given amplitude, assuming no 
further excitation. In that case, the amplitude starts de- 
caying due to strong nonlinear damping, until it reaches 
a threshold where only slow damping occurs. 

In order to detect such an event, one could either 
search in short segments of the detector data for the ini- 
tial strong signal, or search for the tail by investigating 
long time intervals. Since any long tail has an amplitude 
below the onset of strong damping, we already know the 



FIG. 29: Like Figure [28] but for axisymmetric modes. The 
errorbar of the F mode is only correct under the assumption 
that cancellation effects do not increase the error in compar- 
ison to nonradial modes; see Sec. |II C[ 
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FIG. 30: Like Figure |28| but for amplitudes such that the 
damping timescale satisfies r = 20 Tp, where Tp is the oscil- 
lation period of the F mode of the corresponding model. 



maximum strain amplitudes, which are given in Figure [28] 
and |29[ The big unknown is the number of available cy- 
cles, which depends on the residual damping timescale 
of the tail, and thus cannot be computed by nonlinear 
evolution. 

To estimate the detectability of the initial peak, we 
neglect the tail, assume a pure exponential decay, and 
limit the integration time to = 1 s. We define an 
effective strain amplitude by 



^cff = A\ 



1 - 6-27^"/^ 

1 - e-2/(/.^) 



(40) 



where fi is the mode frequency in the inertial frame, and 
T is the initial damping timescale, corrected for the nu- 
merical damping. The enhancement factor in Eq. (40) 
was taken from |52| . For low amplitudes where we can- 
not disentangle numerical and physical damping, we set 
r = 1000 s. 

The results are plotted in Figure [31] As one can see, 
the increase in amplitude after the onset of nonlinearity 
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FIG. 31: Effective GW strain amplitude Aad defined by 
Eq. ( |40[ | versus perturbation amplitude Ai for the nonax- 
isymmetric oscillations. 

competes with the rapid decrease in signal duration, and 
in most cases the latter wins. The effect becomes more 
pronounced when increasing the integration time. When 
considering the tail again, we find that the detectabil- 
ity does not increase significantly for perturbation am- 
pHtudes above the onset of nonlinearity. 



V. SUMMARY AND DISCUSSION 

In our study, we investigated strong nonlinear damping 
effects in high amplitude oscillations of various neutron 
star models, and the limits imposed on the gravitational 
wave signal. For this, we used a relativistic evolution 
code, but kept the spacetime fixed. 

Our neutron star models are rigidly rotating ideal fluid 
configurations with polytropic EOSs of different stiffness, 
and typical neutron star masses in the range 1.4-1.9 so- 
lar masses. The oscillation modes we considered are the 
axisymmetric i^and ^/o modes, the corotating ^/_2 mode, 
and the counter-rotating mode. 

Our main interest is the mode of one model 
(MA65), which can be excited by the CFS instability in 
full GR, and which is therefore a candidate for detectable 
gravitational waves from rapidly rotating neutron stars. 

We found different damping mechanisms for axisym- 
metric and nonaxisymmctric modes. Axisymmetric os- 
cillations of high initial amplitude are rapidly damped 
by formation of shocks in the outer layers of the star, un- 
til the amplitude falls below a certain threshold. Below 
the threshold we observed no further energy dissipation 
on the timescales of our evolutions. 

Although we are using a one-parametric EOS, our re- 
sults accurately determine the onset of shock formation, 
and provide a first estimate for the magnitude of en- 
ergy dissipation by shocks. We find that the damping 
timescale decreases with increasing amplitude, down to 
a few oscillation cycles. The damping timescale and the 
threshold vary by one order of magnitude for the different 



EOSs we investigated. The stiffer the EOS, the higher 
amplitudes are possible before shocks form, and also the 
damping observed at the highest amplitudes is faster for 
stiffer EOSs. 

With increasing rotation rate, nonlinear effects start 
at lower amplitudes. However, the influence of rota- 
tion on the damping of axisymmetric modes is relatively 
weak until close to the Kepler limit, where the damping 
threshold decreases to zero. At the Kepler limit, damp- 
ing by mass shedding occurs. We explain this behavior 
by the fact that the material at the equator is bound less 
strongly for faster rotation. 

For the nonaxisymmctric oscillations, damping is 
mainly caused by wave breaking and nonlinear mode cou- 
pling effects, while shock formation only occurs in low 
density regions at the stellar surface. As for axisymmet- 
ric modes, nonlinear damping sets in at lower amplitudes 
when increasing the rotation rate or decreasing the stiff- 
ness of the EOS. 

Since wave breaking is a surface effect, we feel that an 
accurate description of the damping at high amplitudes 
requires an improved numerical treatment of the surface 
and probably a better physical model as well, considering 
e.g. a hot envelope for a proto-neutron star, a solid crust 
for a cold star (although in that case, one would also have 
to add superfluidity) , and the influence of the magnetic 
field. 

Significant mode coupling was only observed for 
rapidly rotating models, where the main mechanism 
seems to be a 1:2 resonance with high order m = ±2 iner- 
tial modes. For CFS-unstable models, such a resonance 
is always possible. Further, the coupling is most proba- 
bly active already at lower amplitudes, but on timescales 
too long for nonlinear study. Since the growth time of 
CFS-unstable modes is longer as well, it will be necessary 
to investigate the coupling at lower amplitudes perturba- 
tively to determine the saturation amplitudes. Neverthe- 
less, we were able to provide upper limits. 

Next, we computed the gravitational wave strain and 
luminosity of several modes in the linear regime. The 
CFS-unstable mode was found to be a weak emitter; 
at the same amplitude, the strain is thirty times lower 
than for the axisymmetric ^/q mode. This is mainly due 
to the low frequency of the ^/2 mode in the inertial frame. 
To construct models for which this mode is more lumi- 
nous, one needs to further separate the Kepler limit from 
the neutral point (where the ^/2 mode has zero inertial 
frame frequency). This can either be done by increas- 
ing the stiffness of the EOS, or probably by introducing 
differential rotation. Since our CFS-unstable model al- 
ready has a relatively stiff EOS (F = 2.46), we consider 
the latter more realistic. 

We also found that the quasiradial F modes of rapidly 
rotating models emit gravitational waves almost as ef- 
ficient as the ^/o modes. For the CFS-unstable model 
MA65, the F mode emits stronger gravitational waves 
than the mode at the same amplitude. The reason is 
that the eigenfunction of the F mode has a considerable 
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quadrupole moment due to the oblateness of the star, as 
well as the high oscillation frequency. 

Further we investigated the gravitational waves pro- 
duced by oscillations in the nonlinear regime. The max- 
imum strain we observed usually agrees with the values 
extrapolated from the linear regime up to a factor of 3 
(2 for the axisymmetric modes). For the CFS-unstable 
mode however, we found that the signal is dominated 
by secondary modes. We consider those feature 
of our particular setup; their amplitudes most likely de- 
pends on the exact way the main mode was excited. To 
compute the strain for the case of a mode grown to 
high amplitude by the CFS instability, it will be neces- 
sary to model carefully the growth of secondary modes, 
especially at higher frequencies. Secondary modes domi- 
nate the strain only for the CFS-unstable mode, be- 
cause it is the weakest emitter. For models where it has 
a higher frequency, secondary modes are probably less 
important. 

Finally, wc combined the results on strong damping 
and luminosity to establish upper limits on the GW de- 
tectability of various modes. There are two astrophys- 
ically relevant cases: a continuous emission by CFS- 
unstable modes at saturation amplitude, and stable 
modes excited by some violent event. 

From the upper limit for the saturation amplitude of 
the CFS instability, we obtain an upper limit on cor- 
responding gravitational strain. The detcctability of a 
continuous signal depends on the integration time. If we 
consider searches using time windows containing at least 
a few hundred cycles, the ^/2 mode (of model MA65) 
oscillating at the upper limit will be detectable at dis- 
tances around 10 Mpc with advanced LIGO. This dis- 
tance scales approximately quadratically with the mode 
frequency in the inertial frame. Current searches for con- 
tinuous signals in the 5th LIGO science run are not sen- 
sitive to the above source since they exclude the high 
spindown rates caused by gravitational waves of such 
strength. We stress again that the actual saturation am- 
plitude may be considerably lower than the upper limits 



we provided. We also note that we did not take into 
account the secondary excited modes; a careful study of 
the couplings during the growth phase might reveal an 
enhanced detcctability. 

For stable modes excited by a sudden violent event, we 
find that the detcctability does not significantly increase 
for amplitudes above the onset of strong damping, since 
the oscillation is quickly damped below the threshold and 
stays there. The detectability thus depends on the efi^ec- 
tive number of wave cycles in the tail, and thus on the 
residual slow damping at amplitudes near the thresh- 
old. If this information becomes available for a given 
mode, e.g. using perturbative techniques, the detectabil- 
ity can be computed from the luminosities we provided. 
However, we find that the ^/g modes oscillating near the 
strong damping threshold are detectable with advanced 
LIGO at least up to distances of 10 Mpc. The strong 
damping effects thus impose only weak limits on the de- 
tectability of those; for most processes, the important 
factor will probably be the strength of the excitation. 

It is worth mentioning that although some of the ques- 
tions we rise can only be computed with perturbation 
theory, only a nonlinear study like ours can provide the 
necessary limits for the amplitudes at which perturbation 
theory can be safely applied. Given that the physically 
most interesting amplitudes arc at the transition from 
linear to the nonlinear regime, we feel that it is neces- 
sary to combine linear and nonlinear methods to obtain 
accurate results on the detectability of oscillation modes. 
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